home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / MATHPACK.PAS < prev    next >
Pascal/Delphi Source File  |  1985-05-31  |  13KB  |  453 lines
  1. program mathpak;
  2. {(C) Australian Computer Society 1984}
  3. {based on paper A Compact Mathematical Function Package
  4.   by A.P. Clarke and W. Marwood,
  5.   The Australian Computer Journal,
  6.   Vol 16, No 3, August 1984, pp 107 to 114
  7.   entered by W F McGee 1 Dec 1984}
  8. const e=1.0e4;
  9.       sqrt_pi=1.772453850905;
  10.       sqrt2  =1.414213562373;
  11.       eul    =0.577215664901;
  12.       sqrt12 =3.464101615137;
  13.       ln2    =0.693147180559;
  14.       ln10   =2.302585092994;
  15.       ln17   =2.833213344056;
  16.       TERMMAX=500;
  17. type available=(sif,cif,e1f,jnf,mbif,lagf,qf,phif,erff,erfcf,chgf2f,sxf,cxf,qmnf,
  18.                 pmnf,hf,gmf,hgff,chgf1f,ghgff,cf,sf,expf,lf,l10f,done);
  19.      string25=string[25];
  20. var prec:real;
  21.     numparm:array[available]of integer;
  22.     name:array[available]of string25;
  23.     sig_dig:integer;
  24.     sig_dig_flag:boolean;
  25.     sig_dig_real,Ubound,lngmb,gmb,mc_dig:real;
  26.     x:real;
  27. {**********************************************************}
  28. function lnn(x:real):real;
  29. {natural log with zero exception}
  30. begin
  31.   if x=0 then lnn:=-86 else lnn:=ln(x)
  32. {IMPLEMENTATION DEPENDENT  38*ln10}
  33. end;
  34. {***********************************************************}
  35. procedure sum_error(var1,dig1,var2,dig2:real);
  36. var error:real;
  37. begin
  38.   error:=abs(var1)*(exp(-dig1*ln10)+prec) +
  39.          abs(var2)*(exp(-dig2*ln10)+prec);
  40.   sig_dig:=trunc(ln(abs((var1+var2)/error))/ln10)
  41. end;
  42. {***********************************************************}
  43. function ghgf(p,q:integer;a1,a2,a3,a4,c1,c2,c3,c4,x:real):real;
  44. const pqmax=4;
  45. type series=array[0..TERMMAX] of real;
  46.      poch=array[1..pqmax]of real;
  47. var a,c:poch;
  48.     sum,y,abs_y,abs_last_y,sigma:real;
  49.     pqm,pq,k,i,j:integer;
  50.     convergent:boolean;
  51.     ser:series;
  52.  
  53.   procedure accuracy;
  54.     function sig(terms:integer):real;
  55.     var i:integer;psum,s1,s2,weight:real;
  56.     begin
  57.       psum:=0; s1:=sum; s2:=ser[1];weight:=2*(p+q+1);
  58.       for i:=1 to terms-1 do begin
  59.         s1:=s1-ser[i-1];s2:=s2+ser[i+1];
  60.         psum:=psum+sqr(s1)+weight*sqr(s2)
  61.       end;
  62.       sig:=sqrt(psum+sqr(s1)+weight*sqr(ser[terms]))
  63.     end;
  64.  
  65.   begin
  66.     sigma:=prec/sqrt12*sig(k-1);
  67.     if sigma=0 then sig_dig_real:=mc_dig
  68.     else if sum=0 then sig_dig_real:=0
  69.     else sig_dig_real:=ln(abs(sum/sigma))/ln10;
  70.     if sig_dig_real<=0.0 then sig_dig_real:=0.0;
  71.     if sig_dig_real> mc_dig then sig_dig_real:=mc_dig;
  72.     sig_dig:=trunc(sig_dig_real)
  73.   end;
  74. begin
  75.   a[1]:=a1; a[2]:=a2; a[3]:=a3; a[4]:=a4;
  76.   c[1]:=c1; c[2]:=c2; c[3]:=c3; c[4]:=c4;
  77.   if p>q then pqm:=p else pqm:=q;
  78.   y:=1.0; sum:=1.0; k:=1; ser[0]:=1.0; pq:=p+q;
  79.   convergent:=false; abs_last_y:=1.0;
  80.   repeat
  81.     for i:=1 to pqm do y:=y*a[i]/c[i];
  82.     for i:=1 to p do a[i]:=a[i]+1;
  83.     for i:=1 to q do c[i]:=c[i]+1;
  84.     y:=y*x/k; abs_y:=abs(y);
  85.     if (abs_y < abs_last_y) then convergent:=true;
  86.     if (abs_y > abs_last_y) and convergent then  y:=0;
  87.     sum:=sum+y; abs_last_y:=abs(y);ser[k]:=y; k:=k+1
  88.   until (sum=sum-y) or (k>499);
  89.   if (k>499) then writeln('not convergent after 499 iterations');
  90.   if sig_dig_flag then accuracy;
  91.   ghgf:=sum
  92. end;
  93.  
  94. procedure init_consts;
  95. begin
  96.   sig_dig:=0;Ubound:=10; gmb:=17; lngmb:=ln17;
  97.   mc_dig:=ln(1/prec)/ln10
  98. end;
  99.  
  100. function chgf1(a,c,x:real):real;
  101. {confluent hypergeometric function of the first kind}
  102. begin
  103.   if x<0.0 then chgf1:=exp(x)*chgf1(c-a,c,-x)
  104.   else chgf1:=ghgf(1,1,a,1,1,1,c,1,1,1,x)
  105. end;
  106.  
  107. function si(x:real):real;
  108. {sine integral}
  109. begin
  110.   si:=x*ghgf(1,2,0.5,1,1,1,1.5,1.5,1,1,-x*x/4.0)
  111. end;
  112.  
  113. function e1(x:real):real;
  114. {Exponential integral}
  115. var var1,var2:real;
  116. begin
  117.   var1:=x*ghgf(2,2,1,1,1,1,2,2,1,1,-x); var2:=-ln(x);
  118.   e1:=-eul+var1+var2;
  119.   if sig_dig_flag then begin
  120.     sum_error(-eul,mc_dig,var1,sig_dig_real);
  121.     sum_error(-eul+var1,sig_dig_real,var2,mc_dig)
  122.   end
  123. end;
  124.  
  125. function ci(x:real):real;
  126. {Cosine integral}
  127. var var1,var2:real;
  128. begin
  129.   var1:=-sqr(x/2)*ghgf(2,3,1,1,1,1,1.5,2,2,1,-sqr(x/2)); var2:=ln(x);
  130.   ci:=eul+var1+var2;
  131.   if sig_dig_flag then begin
  132.     sum_error(eul,mc_dig,var1,sig_dig_real);
  133.     sum_error(eul+var1,sig_dig_real,var2,mc_dig)
  134.   end
  135. end;
  136.  
  137. function hgf(a1,a2,c,x:real):real;
  138. {Gauss hypergeometric function}
  139. begin
  140.   hgf:=ghgf(2,1,a1,a2,1,1,c,1,1,1,x)
  141. end;
  142.  
  143. function gm(x:real):real;
  144. {Gamma function}
  145. var temp_flag:boolean;
  146. begin
  147.   temp_flag:=sig_dig_flag; sig_dig_flag:=false;
  148.   if x=0.0 then gm:=1.0
  149.   else gm:=exp(x*lngmb)/x*chgf1(x,x+1,-gmb)+
  150.     exp((x-1)*lngmb-gmb)*ghgf(2,0,1-x,1,1,1,1,1,1,1,-1.0/gmb);
  151.   sig_dig_flag:=temp_flag
  152. end;
  153.  
  154. function h(n:integer;x:real):real;
  155. {Hermite polynomials}
  156. var sign:integer;
  157. begin
  158.   if odd(n div 2) then sign:=-1 else sign:=1;
  159.   if (n mod 2)=0 then h:=sign*gm(n+1)/gm(n/2+1)*chgf1(-n/2,0.5,x*x)
  160.   else h:=sign*gm(n+1)/gm((n+1)/2)*2*x*chgf1(-(n-1)/2,1.5,x*x)
  161. end;
  162.  
  163. function pmn(m,n:integer;x:real):real;
  164. {Associated Legendre polynomial of the first kind}
  165. var sign:real;
  166. begin
  167.   if abs(x)>1 then begin
  168.     writeln('usage: pmn needs abs(argument)<1');
  169.     pmn:=0.0
  170.   end
  171.   else begin
  172.     if odd(m) then sign:=-1.0 else sign:=1.0;
  173.     pmn:=gm(m+n+1)/(gm(m+1)*gm(n-m+1))*exp(ln(1-sqr(x))*m/2-ln2*m)*
  174.          sign*ghgf(2,1,m-n,m+n+1,1,1,1+m,1,1,1,(1-x)/2)
  175.   end
  176. end;
  177.  
  178. function qmn(m,n:integer; x:real):real;
  179. {Associated Legendre polynomial of the second kind}
  180. var sign:real;
  181.     temp1:real;
  182. begin
  183.   if abs(x)<1 then begin
  184.     writeln('usage:qmn needs abs(argument) >1');
  185.     qmn:=0.0
  186.   end
  187.   else begin
  188.     if odd(m) then sign:=-1.0 else sign:=1.0;
  189.     temp1:=(lnn(x*x-1.0)*m/2.0)-lnn(x)*(n+m+1.0);
  190.     qmn:=exp(temp1-ln2*(n+1))*sqrt_pi*
  191.       gm(n+m+1)/gm(n+1.5)*sign*
  192.       ghgf(2,1,(1+n+m)/2,(2+m+n)/2,1,1,n+1.5,1,1,1,1/(x*x))
  193.   end
  194. end;
  195.  
  196. function cx(x:real):real;
  197. {Fresnel integral of the first kind}
  198. begin
  199.   cx:=x*ghgf(1,2,0.25,1,1,1,0.5,1.25,1,1,-sqr(pi*sqr(x/2)))
  200. end;
  201.  
  202. function sx(x:real):real;
  203. {Fresnel integral of the second kind}
  204. begin
  205.   sx:=pi*x*x*x/6*ghgf(1,2,0.75,1,1,1,1.5,1.75,1,1,-sqr(pi*sqr(x/2)))
  206. end;
  207.  
  208. function chgf2(a,b,x:real):real;
  209. {confluent hypergeometric function of the second kind}
  210. var temp1,temp2,sig1:real;
  211. begin
  212.   if abs(x)<Ubound then begin
  213.     temp1:=chgf1(a,b,x)/gm(1+a-b)/gm(b);sig1:=sig_dig_real;
  214.     temp2:=chgf1(1+a-b,2-b,x)*exp((1-b)*lnn(x))/gm(a)/gm(2-b);
  215.     chgf2:=(temp1-temp2)*pi/sin(pi*b);
  216.     if sig_dig_flag then sum_error(temp1,sig1,-temp2,sig_dig_real)
  217.   end
  218.   else chgf2:=exp(-a*lnn(x))*ghgf(2,0,a,1+a-b,1,1,1,1,1,1,-1.0/x)
  219. {Asymptotic representation for large abs(x)}
  220. end;
  221.  
  222. function erfc(x:real):real;
  223. {Complementary error function}
  224. begin
  225.  if x<0 then erfc:=2.0-exp(-x*x)*chgf2(0.5,0.5,x*x)/sqrt_pi
  226.  else erfc:=exp(-x*x)*chgf2(0.5,0.5,x*x)/sqrt_pi
  227. end;
  228.  
  229. function erf(z:real):real;
  230. {Error function}
  231. begin
  232.   if abs(z)<3.0 then erf:=z*chgf1(0.5,1.5,-z*z)*2/sqrt_pi
  233.   else erf:=1.0-erfc(z)
  234. end;
  235.  
  236. function phi(z:real):real;
  237. {Normal probability integral}
  238. var var1:real;
  239. begin
  240.   var1:=erf(z/sqrt2);phi:=0.5+0.5*var1;
  241.   if sig_dig_flag then sum_error(0.5,mc_dig,0.5*var1,sig_dig_real)
  242. end;
  243.  
  244. function q(x:real):real;
  245. {Complementary Normal probability integral}
  246. begin
  247.   q:=0.5*erfc(x/sqrt2)
  248. end;
  249.  
  250. function lag(alpha,n,z:real):real;
  251. {Laguerre polynomials}
  252. begin
  253.   lag:=gm(alpha+n+1)/(gm(n+1)*gm(alpha+1))*chgf1(-n,alpha+1,z)
  254. end;
  255.  
  256. function mbi(nn,x:real):real;
  257. {Modified Bessel function}
  258. var n:real;
  259. begin
  260.   n:=nn;
  261.   if round(n)=n then n:=abs(n);
  262.   mbi:=exp(n*lnn(x/2))/gm(n+1)*ghgf(0,1,1,1,1,1,1+n,1,1,1,sqr(x/2))
  263. end;
  264.  
  265. function jn(n,x:real):real;
  266. {Bessel function}
  267. begin
  268.   jn:=exp(n*lnn(x/2))/gm(n+1)*ghgf(0,1,1,1,1,1,1+n,1,1,1,-sqr(x/2))
  269. end;
  270.  
  271. procedure initfunctions;
  272. begin
  273.   name[sif]:='Sine-Integral';
  274.   name[cif]:='Cosine-Integral';
  275.   name[e1f]:='Exponential-Integral';
  276.   name[jnf]:='Bessel';
  277.   name[mbif]:='Modified-Bessel';
  278.   name[lagf]:='Laguerre';
  279.   name[qf]:='Complementary-Normal';
  280.   name[phif]:='Normal-Probability';
  281.   name[erff]:='Error';
  282.   name[erfcf]:='Complementary-Error';
  283.   name[chgf2f]:='Confluent-2nd-type';
  284.   name[sxf]:='2nd-Fresnel-Integral';
  285.   name[cxf]:='1st-Fresnel-Integral';
  286.   name[qmnf]:='2nd-Type-Legendre-(Assc)';
  287.   name[pmnf]:='1st-Type-Legendre-(Assc)';
  288.   name[hf]:='Hermite';
  289.   name[gmf]:='Gamma';
  290.   name[hgff]:='Gauss-Hypergeometric';
  291.   name[chgf1f]:='Confluent-1st-Type';
  292.   name[ghgff]:='Generalized-Hyper';
  293.   name[cf]:='Cosine';
  294.   name[sf]:='Sine';
  295.   name[lf]:='Ln';
  296.   name[l10f]:='Log10';
  297.   name[expf]:='Exp';
  298.   name[done]:='DONE';
  299.  
  300.   numparm[sif]:=1;
  301.   numparm[cif]:=1;
  302.   numparm[e1f]:=1;
  303.   numparm[jnf]:=2;
  304.   numparm[mbif]:=2;
  305.   numparm[lagf]:=3;
  306.   numparm[qf]:=1;
  307.   numparm[phif]:=1;
  308.   numparm[erff]:=1;
  309.   numparm[erfcf]:=1;
  310.   numparm[chgf2f]:=3;
  311.   numparm[sxf]:=1;
  312.   numparm[cxf]:=1;
  313.   numparm[qmnf]:=3;
  314.   numparm[pmnf]:=3;
  315.   numparm[hf]:=2;
  316.   numparm[gmf]:=1;
  317.   numparm[hgff]:=4;
  318.   numparm[chgf1f]:=3;
  319.   numparm[ghgff]:=11;
  320.   numparm[cf]:=1;
  321.   numparm[sf]:=1;
  322.   numparm[lf]:=1;
  323.   numparm[l10f]:=1;
  324.   numparm[expf]:=1;
  325.   numparm[done]:=0;
  326. end;
  327.  
  328. function calculate:available;
  329. var i:integer;r:real;index:available;
  330.     getstring:string25;
  331.     vari:array[1..11]of real;
  332. procedure getvar(num:integer);
  333. var i:integer;
  334. begin
  335.   writeln('Entry of parameters');
  336.   if num=0 then
  337.   else for i:=1 to num do begin
  338.     write('Parameter [ ',i,'] is ');readln(vari[i])
  339.   end
  340. end;
  341. begin
  342.   index:=sif;
  343.   repeat
  344.     write(name[index]:30);
  345.     index:=succ(index);
  346.     if index<>done then begin
  347.       writeln(name[index]:30);index:=succ(index)
  348.     end
  349.   until index=done;
  350.   writeln;writeln('Enter function desired');
  351.   readln(getstring);
  352.   index:=sif;
  353.   while not ((index=done)or(name[index]=getstring)) do index:=succ(index);
  354.   getvar(numparm[index]);
  355.   writeln;writeln('Working....');
  356.   if index=sif then r:=si(vari[1]);
  357.   if index=cif then r:=ci(vari[1]);
  358.   if index=e1f then r:=e1(vari[1]);
  359.   if index=jnf then r:=jn(round(vari[1]),vari[2]);
  360.   if index=mbif then r:=mbi(round(vari[1]),vari[2]);
  361.   if index=lagf then r:=lag(vari[1],vari[2],vari[3]);
  362.   if index=qf then r:=q(vari[1]);
  363.   if index=phif then r:=phi(vari[1]);
  364.   if index=erff then r:=erf(vari[1]);
  365.   if index=erfcf then r:=erfc(vari[1]);
  366.   if index=chgf2f then r:=chgf2(vari[1],vari[2],vari[3]);
  367.   if index=sxf then r:=sx(vari[1]);
  368.   if index=cxf then r:=cx(vari[1]);
  369.   if index=qmnf then r:=qmn(round(vari[1]),round(vari[2]),vari[3]);
  370.   if index=pmnf then r:=pmn(round(vari[1]),round(vari[2]),vari[3]);
  371.   if index=hf then r:=h(round(vari[1]),vari[2]);
  372.   if index=gmf then r:=gm(vari[1]);
  373.   if index=hgff then r:=hgf(vari[1],vari[2],vari[3],vari[4]);
  374.   if index=chgf1f then r:=chgf1(vari[1],vari[2],vari[3]);
  375.   if index=ghgff then r:=ghgf(round(vari[1]),round(vari[2]),vari[3],vari[4],
  376.                           vari[5],vari[6],vari[7],vari[8],vari[9],
  377.                           vari[10],vari[11]);
  378.   if index=cf then r:=cos(vari[1]);
  379.   if index=sf then r:=sin(vari[1]);
  380.   if index=lf then r:=ln(vari[1]);
  381.   if index=expf then r:=exp(vari[1]);
  382.   if index=l10f then r:=ln(vari[1])/ln10;
  383.   if index=done then r:=0.0;
  384.   clrscr;
  385.   write(name[index],'[');
  386.   if numparm[index]<>0 then begin
  387.     for i:=1 to numparm[index] do begin
  388.      write(vari[i]:12:6);if i<>numparm[index] then write(' , ')
  389.    end
  390.   end;
  391.   write(' ]= ');writeln(r);
  392.   writeln;writeln('significant digits =  ',sig_dig);
  393.   calculate:=index
  394. end;
  395.  
  396. function precision:real;
  397. var x,y:real;
  398. begin
  399.   y:=1.0;
  400.   x:=1.0;
  401.   repeat
  402.     x:=x+y;
  403.     y:=10.0*y
  404.   until x+1.0=x;
  405.   x:=y/100.0;
  406.   repeat
  407.     x:=x+y/1000.0
  408.   until x+1.0=x;
  409.   x:=x-y/1000;
  410.   precision:=1.0/x
  411. end;
  412. begin
  413.   clrscr;
  414.   writeln('FUNCTION PACKAGE BY CLARKE AND MARWOOD');
  415.   writeln('for TURBO PASCAL  peak exponent 10E(+/- 38)');
  416.   writeln('ignore digits precision for sine,cos,exp,ln,log10');
  417.   prec:=precision;
  418.   writeln('PRECISION =',prec);
  419.   init_consts;
  420.   sig_dig_flag:=true;
  421.   initfunctions;
  422.   repeat
  423.   until calculate=done
  424. end.
  425. function Q(x:real):real;
  426. var temp,factor,con,z,q1:real;
  427. begin
  428.   con:=2.506628274;
  429.   if abs(x)>12.5 then q1:=0.0
  430.   else begin
  431.     temp:=sqr(x);
  432.     z:=exp(-temp/2.0);
  433.     factor:=abs(x/2)+sqrt(1.0+temp/4.0);
  434.     q1:=z/(con*factor+(2.0-con)*z)
  435.   end;
  436.   if x>=0.0 then Q:=q1 else Q:=1.0-q1;
  437. end;
  438.  
  439. function arcq(qd:real):real;
  440. var x0,x:real;
  441. begin
  442.   if abs(qd)>=1.0 then writeln('out of range')
  443.   else begin
  444.     if qd>0.5 then x0:=0 else x0:=sqrt(-2*ln(qd))-1.0;
  445.     repeat
  446.       x:=x0;
  447.       x0:=x0-(q(x0)-qd)/(-exp(-sqr(x0)/2)/sqrt(2*pi));
  448.     until abs(x0-x)<0.00000001
  449.   end;
  450.   arcq:=x0
  451. end;
  452.      x:=x0;
  453.       x0:=x0-(q(x0)-qd)/